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Abstract 

We perform density functional molecular dynamics simulations of liquid and 
solid MgSi0 3 in the pressure range of 120—1600 GPa and for temperatures up 
to 20000 K in order to provide new insight into the nature of the liquid-liquid 
phase transition that was recently predicted on the basis of decaying laser 
shock wave experiments [Phys. Rev. Lett. 108 (2012) 065701]. However, 
our simulations did not show any signature of a phase transition in the liquid 
phase. We derive the equation of state for the liquid and solid phases and 
compute the shock Hugoniot curves. We discuss different thermodynamic 
functions and by explore alternative interpretations of the experimental find- 
ings. 

Keywords: density functional molecular dynamics, ab initio simulations, 
high pressure, shock wave experiments 



1. Introduction 

The recent work by D. K. Spaulding et al. [l[ reported results from decay- 
ing laser shock wave experiments which provided evidence of a liquid-liquid 
phase transition in MgSi03 at megabar pressures. The authors measured 
a reversal in the shock velocity and thermal emission and interpreted their 
findings in terms of a liquid-liquid phase transition that occurs when the 
sample changes from a high density to a low density fluid state during shock 
decay. Decaying shock experiments are a new experimental technique that 
allows one to map out a whole segment of the Hugoniot curve with a single 
shock wave experiment. 
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First-principles computer simulations have a long tradition of character 



and 




izing materials at extreme pressures 0, 0, 0] and temperatures 0, 0, 
of contributing to the interpretation of shock wave experiments 
The goal of this particular paper is to perform density functional molecular 
dynamics simulations (DFT-MD) of dense liquids jsj in order to verify the 
predictions of a liquid-liquid phase transition by Spaulding et al. First 
order transitions in liquids are unusual but have been seen in experiments 
on phosphorus and in simulations of dense hydrogen 
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2. Simulation Parameters 



All simulations were performed with the VASP code [12] with pseudopo 



tentials of the projector-augmented wave type 13], a cut-off for the expansion 



of the plane wave basis set for the wave functions of 500 eV, and the PBE 



exchange-correlation functional [141 ] . The Brillioun zone was only sampled 
with the T point to allow for extended and efficient MD simulations. Our 
simulations lasted between 1 and 20 ps and used a small time step of 0.2 fs. 



The electronic states were populated according to the Mermin functional [15 

We used an orthorhombic supercell with 60 atoms that we constructed by 
tripling the unit cell of the post-perovskite (PPV) structure that we relaxed 
at 120 GPa. The system was heated and melted using a Nose thermostat. 
We then explored the liquid state by scaling the velocities and changing the 
density accordingly. 

We also performed heat-until-it-melts simulations with 60 atoms starting 
a perfect PPV crystal at hydrostatic conditions. We then gradually increased 
the temperature in a fixed cell geometry. We also performed a large number 
of solid simulations at constant temperature and different densities. 



3. Results and Discussion 

The P-T conditions of our simulations are shown in Fig. Q] and the re- 
sulting equation of state (EOS) for the liquid and the post-perovskite solid 
phases are given in Tables [1] and [2j Initially we focused our work on 7500 K 
and 12500 K in order to maximize the likelihood of directly observing a phase 
transition in the liquid during the MD simulations. However, all the ther- 
modynamic functions, that we analyzed, varied in a perfectly smooth way as 
function of temperature and pressure. Figure [2] shows the volume and the 
internal energy as function of pressure. The zero of energy was set equal to 



2 




Figure 1: In the left diagram, the open diamond and circles indicate the pressure- 
temperature conditions of our solid and liquid MgSiC>3 DFT-MD simulations, respectively. 
The thin black dashed lines show the P-T paths of the solid portions of our heat-until- 
it-melts simulations. The computed solid and liquid shock Hugoniot curves are shown in 
blue and red color corresponding to initial densities of 3.22 and 2.74 gcm~ 3 to match 
those of the single crystal (blue) and glass (red) materials used in the experiments 
The Hugoniot curves are repeated in black and white on the right for comparison with the 
decaying shock measurements From the discontinuities in the experimental Hugoniot 
curves, the phase boundary between a high and low density liquid was inferred (dash- 
dotted line). The computed Hugoniot curve for solid MgSiOs agree well with earlier gas 



gun shock measurements (diamonds) [16l Il7j . Three earlier predictions of the melting line 
are included [H, [TEH^- 
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the energy of a collection of isolated atoms. We conclude from this analysis, 
if there exists a phase transition in the liquid phase then ab initio simulations 
of 60 atoms are unlikely to spontaneously transform into a different phase. 




Figure 2: Volume and internal energy as function of pressure derived from DFT-MD 
simulations of liquid MgSi03. Both functions are smooth and exhibit no indication of a 
phase transition. 



We computed the shock Hugoniot curve from the usual relation 21], 



H = E - E + 1 -{V - V )(P + P ) = 0, (1) 



where the initial volume, Vo=51.77 A 3 , and initial internal energy, Eq 
—35.914 eV per MgSiC>3 formula unit (FU), were taken from a DFT calcula- 
tion of enstatite at 3.22 gem -3 that we performed. We used the same E to 
derive the Hugoniot curves for a glass sample with an initial density of 2.74 
gem" 3 . Pq was assumed to be zero because Pq <C P is satisfied for all final 
states. 

The resulting DFT-MD Hugoniot curves for liquid and crystalline MgSi03 
are shown in Fig. [TJ Our results for the solid phase agree well with the 



earlier shock measurements using a gas gun [16|, [17j . However, the agreement 
with the decaying laser shock results is not satisfactory. A portion of our 
computed Hugoniot curve for the solid at po = 3.22 gem" 3 overlaps with the 
measurements in the temperature range from 7000 to 9000 K. Also at 450 
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Table 1: Temperature, MD simulation time, density, pressure and internal energy from 
our DFT— MD simulations of liquid MgSiOg. The la uncertainties of the trailing digits 
are given in brackets. 



T(K) 


time (ps) 


p (gem" 3 ) 


P (GPa) 


E/FXJ (eV) 
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5.00 
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Table 2: Temperature, MD simulation time, density, pressure and internal energy from 
our DFT-MD simulations of crystalline post-perovskite MgSiC>3. The la uncertainties of 
the trailing digits are given in brackets. 
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GPa and 12500 K, the calculated Hugoniot curve for the liquid agrees with 
the measurements but the slope of the theoretical Hugoniot curve is different 
and it shows no sign of a phase transition. The agreement for an initial 
density of 2.74 gem -3 is again not favorable. The theoretical Hugoniot the 
liquid passes right through the T-P conditions where a phase transformation 
is predicted based on the experimental data but DFT-MD simulations do 
not exhibit one. 
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Figure 3: Diffusion coefficients of Mg, Si, and O ions as function of pressure for 7500 
and 12500 K inferred from DFT-MD simulations. The lines show fits to an exponential 
function of pressure for each species. 

Next we investigated the possibility of the existence of a superionic phase 
between the solid and the liquid phases where some ions remain stationary 



while others diffuse throughout the material like in a liquid [22J. In Fig. [21 we 
plot the diffusion coefficients as a function of pressure. At 12500 K, all ions 
diffuse at the same rate approximately. At 7500 K, the diffusion is slower as 
one would expect in the fluid near the melting line. The oxygen ions diffuse 
a bit faster than the magnesium and silicon ions but there is no evidence for 
a superionic state. 

We also performed heat-until- it-melts simulations of solid samples (Fig. [1]) 
that we isochorically heated at a rate of 1000 K per picosecond. This ap- 



proach was used to predict the superionic state of water 22] but our sample 
always went into a fully fluid state after a substantial amount of superheat- 
ing. Figure H] shows that the liquid always exhibits a higher pressure than the 
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Figure 4: The pressure-temperature behavior during heat-until-it-melts simulations (lines) 
is compared with isothermal liquid simulations (symbols). The small arrows indicate the 
beginning of the melting transformation in the MD simulation. Liquid samples are found 
to always exhibit a higher pressure than solid ones at the same density. 



solid when compared for the same density and temperature. This means, for 
given temperature and pressure, the solid is always denser than the liquid. 
Therefore the melting line of PPV MgSiOs is expected to have a positive 
Clapeyron slope [23 . 

We performed our simulations with only 60 atoms and all results may 
need to be confirmed with larger simulations that possibly also use more 
k-points. Our longest simulations ran for 20 ps, which is fairly long for 
ab initio calculations. If a new liquid phase emerges, one would expect to 
observe a spontaneous phase transformation, as long as it does not involve 
any large-scale rearrangements of atoms. If liquid MgSiOa phase-separates 
into a MgO-dominated and a SiCVrich fluid, as was predicted for the solid 
at approximately 1000 GPa [24j , we would most likely not be able to observe 
this process in our simulations. One would instead need to perform more 
expensive Gibbs free energy calculations [25|, |26|, \4, |27|, |28j but such a complex 
effort goes beyond this investigation. 

Furthermore, if there existed a new, unknown phase in the MgSiOa phase 
diagram, it does not have to be a liquid. One could imagine a new stable 
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solid phase that introduces a solid-solid-liquid triple point and may lead to an 
increase in the slope of the melting line. At zero temperature, existence of a 
post-post-perovskite phase has been intensely studied with ab initio random 
structure search algorithms and no such phase has been found but a new, 
entropy-stabilized phase may still exist at high temperature. In principle, this 
solid-solid-liquid triple point could also be between perovskite (PV), PPV, 
and liquid phases. This point has not yet been determined neither with 
experimental nor with theoretical means but PV-to-PPV transition pressure 
is known to increases with temperature. 

One may also ask whether there exists an alternative interpretation for 
the experimental observations. Our first recommendation would be to repeat 
those measurements with steady shock wave experiments in order to verify 
the discontinuities on the principal Hugoniot curve of MgSiOa. This may 
require a series of shock experiments with relatively small error bars but it 
would be an important confirmation. 

It is also possible that observed shock velocity reversal is related to the 
melting transition of MgSi03. Figure [T] shows that most of the experimental 
results fall in between the computed Hugoniot curves for the solid and the 
liquid. However, the temperatures are too high for a solid phase to be ther- 
modynamically stable. Phase transition temperatures of 10 000 and 16 000 K 
have been predicted based on the measurements of single- crystal and glass 
material [if, respectively. Nevertheless it is interesting to discuss the behav- 
ior of decaying shock experiments in the presence of a melting transition. 
Some material behind the shock front may freeze during the shock decay. 
This could introduce a secondary wave and thereby affect the behavior of 
the shock front. 

At the beginning of a decaying shock experiments, the sample material is 
compressed to a state of high pressure and high temperature on the principal 
Hugoniot curve. As the shock decays, new material does not reach as high 
pressures and temperatures but stays on the Hugoniot curve, which allows 
one to map out the whole Hugoniot curve with just one shock measurement. 
The question is what happens to the material that had compressed to a higher 
P-T state earlier. One may assume that the whole region behind the shock 
front will equilibrate to a new pressure, which we labeled P* in Fig |5j Any 
parcel of material that was shocked to a high P-T state will adiabatically 



expand 21] to reach P* and slow down to travel at the new and reduced 



particle velocity. Since this expansion is gradual and not associated with 



any shock, one typically also assumes that this expansion is isentropic [21 
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Figure 5: Temperature-pressure diagram with melting lines and shock Hugoniot curves. 
The arrows indicate the adiabatic and isentropic expansion of material during decaying 
shocks. The graph shows a material with normal melting behavior with a positive Clapey- 
ron slope, dT m /dP > 0. 

(arrows in Fig. [5]). While one assumes a new pressure, P*, is established, 
the entire sample behind the shock front is not expected to reach thermal 
equilibrium during the experiment. This leads to the situation where hot 
material at a lower density is pushing on colder material at a higher density. 
This could lead to a fluid dynamic instability of Rayleigh- Taylor type if the 
density contrast becomes too high and the time scale is long enough for such 
an instability to develop. 

It is important to note that Hugoniot curves are steeper in P-T space 
than isentropes, |^| Hug > The shock front on the Hugoniot curve 

would always enter any new thermodynamic phase, that may exist at lower 
temperature, before the hotter material behind the shock front reaches it. 

However, the melting transition may have an effect on the shock propa- 
gation at much higher temperature. Zeldovich and Raizer [2l| showed that 
a shock will split into two waves if the shock Hugoniot curves intersects with 
the Rayleigh line that linearly connects the initial and final states in a P-V 
diagram. If the melting transition would introduce such a splitting, it would 
affect the shock propagation at much higher temperatures and pressures. We 
plotted our liquid and solid Hugoniot curves the P-V and u s -u p diagrams in 
Figures |6] and [7J For each of the Hugoniot curves, let us assume that ma- 
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Figure 6: P-V diagram of the solid and liquid segments of the shock Hugoniot curves for 
two initial densities. The labels indicate the calculated shock temperatures. 



terial melts at a temperature, T m , somewhere in the range spanned by the 



calculation in Refs. [18|, |19|, |20j (see Fig. [T)). If there occurred a sudden and 
complete phase change from liquid to solid during a decaying shock experi- 
ments at T m , it would lead to a abrupt drop in pressure, density, shock and 
particle velocities because intermediate values cannot be realized by either 
phase while satisfying Eq. (Op). 

From Fig. El we conclude that the Hugoniot curves do not intersect with 
the Rayleigh line but the material must assume a mixed state because the 
particle velocity decreases gradually in a decaying shock wave experiment. 
Following the arguments that Greeff et al. [29| applied to the a-to-u; transi- 
tion in solid titanium, the properties of a solid-liquid MgSiOs mixture with 
the solid fraction, A, can be derived from 

E = (l-X)E l (P,T) + XE s (P,T) (2) 
V = (1 - XMP,T) + \V S (P,T) (3) 

for a given pair of T and P on the melting line. Under the assumption of 
thermodynamic equilibrium, A can derived from Eq. ([T]) for given pressure 
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Figure 7: Shock velocity versus particle velocity for the Hugoniot curves from Fig. [6] 

or particle velocity. For high shock velocities, the final state of the material 
is fully molten. As the shock decays, the solid fraction gradually increases 
until the shock is too weak to introduce any melting. 

We now want to address the question what happens as the parcels of solid- 
liquid mixtures expand adiabatically. The solid fraction, A, must change to 
keep the entropy of the mixture, 

S = (l-X)S l (P,T) + XS s (P,T) , (4) 

constant as the pressure decreases. T adjusts to the melting temperature for 
given pressure. If the solid fraction, A, gradually increases during the ex- 
pansion, the particle velocity directly behind the shock front would decrease. 
This is not unusual but the opposite case is much more interesting. If the 
liquid fraction increases during the shock decay, IsL > 0, it would imply 

Or I 1 m 

an increase in particle velocity. As a result the shock wave would split in 
two waves, which would make the analysis of the VISAR signals much more 
difficult. 

We can determine the thermodynamic parameters for such a shock wave 
splitting. During the expansion, the total entropy stays constant, ff L = 0. 



12 



Eq. PJ implies, 



5A_ 



(5, - 5 S ) = (1 - A) 



as* 



+ A 



dS s 



8P rp 

1 f 



(5) 



Since the entropy of the liquid is always greater than that of the solid, a 
positive right term in Eq. (jSj) would imply that the shock wave splits into 
two separate waves as the decay shock experiment enter the regime of solid- 
liquid coexistence. Without a detailed calculations of the entropies in both 
phases, it is not possible to determine whether this occurs in MgSiOs because 
the Clapeyron slope enters into the calculation. However, the entropy is 
accessible with thermodynamic integration techniques and equation of state 
models [20]. So our hypothesis of shock wave splitting can be tested with 
experimental and theoretical techniques. 

Summarizing we can say that whether a Rayleigh- Taylor instability exists, 
a new solid phase appears, the shock wave splits into two waves, or indeed 
a liquid-liquid phase transition is the reason for the observed shock velocity 
reversal remains to be determined with future experiments. Concluding our 
theoretical investigation, we can report that our DFT-MD simulations to 
not support the hypothesis of a liquid-liquid phase transition in MgSiOs. We 
were not able to give an obvious alternative interpretation of the experimen- 
tal findings but we do, however, suggest the measurements be repeated with 
steady shock waves. The experimental findings are very interesting never- 
theless and will lead to a better understanding of magnesiosilicates at high 
temperature. 
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